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The three-dimensional g-state Potts model, forced into co- 
existence by fixing the density of one state, is studied for 
q = 2, 3, 4, and 6. As a function of temperature and number 
of states, we studied the resulting equilibrium droplet shapes. 
A theoretical discussion is given of the interface properties 
at large values of q. We found a roughening transition for 
each of the numbers of states we studied, at temperatures 
that decrease with increasing q, but increase when measured 
as a fraction of the melting temperature. We also found equi- 
librium shapes closely approaching a sphere near the melting 
point, even though the three-dimensional Potts model with 
three or more states does not have a phase transition with a 
diverging length scale at the melting point. 



I. INTRODUCTION 

If a binary mixture of fixed composition is brought 
into a coexistence state it will phase-separate into two 
pure phases, separated by an interface with a shape that 
minimizes the free energy of the system The phase 
with the smaller volume will typically organize itself into 
a compact shape, known as the equilibrium shape. If 
the surface tension (excess free energy per unit inter- 
face area) is isotropic, this equilibrium shape will be a 
sphere. If either of the co-existing phases is crystalline, 
the anisotropic surface tension will lead to an aspheric 
equilibrium crystal shape (ECS). The ECS and the ori- 
entation dependence of the surface tension are intimately 
related. Once the surface tension is known for all orien- 
tations, the so-called Wulff construction |2|,|| allows the 
generation of the ECS. The reverse is also possible, i.e., 
determining the surface tension as a function of orienta- 
tion from the ECS; this is a procedure followed in this 
manuscript. 

The ECS has been the subject of experimental stud- 
ies jiffy]. One striking feature the ECS may show 
is a roughening transition: the disappearance of facets 
(macroscopically flat surfaces) under rising temperature. 
In practice equilibrium shapes are rarely seen and they 
are hard to produce experimentally, but this does not 
mean they are without practical importance. For exam- 
ple the presence or absence of facets in the ECS influ- 
ences growth properties, such as the speed of growth and 
the growth mode, even though growth shapes may dif- 
fer strongly from the ECS. In general, the determination 
of the orientation dependence of the surface tension is 
a tough problem. Experimentally it is very difficult to 



measure, mainly due to equilibration problems. Numer- 
ical studies suffer from the same equilibration problems, 
and in the most common simulation techniques each sur- 
face orientation requires a separate simulation. Theoret- 
ical results only have been obtained for simplified models 
like the BCSOS model 

The prototypical model in which such properties are 
studied numerically, is the conserved-order-parameter 
Ising model with nearest-neighbor interactions, in which 
the total magnetization is kept constant. In the two- 
dimensional Ising model, the ECS is a square at zero tem- 
perature, but at any finite temperature the ECS looses all 
flat faces. It gradually changes into a circle when the tem- 
perature is approaching the critical temperature. The be- 
havior is richer in the three-dimensional Ising model. At 
zero temperature, the ECS is a cube. At finite temper- 
atures below the roughening temperature Tr, the ECS 
still has macroscopically flat faces but the corners as well 
as the edges are rounded. Above the roughening temper- 
ature, the ECS does no longer feature macroscopically 
flat faces. If the temperature is increased from the rough- 
ening temperature to the critical temperature, the ECS 
gradually approaches a sphere. 

An extension of the Ising model is the Potts model, 
defined by the Hamiltonian 
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in which J is the coupling constant, S denotes the Kro- 
necker delta function, and the summation runs over all 
pairs of nearest-neighbor sites, each having a spin with 
value a = 1 . . . q. Note that the two-state Potts model 
is equivalent to the Ising model. The topic of this 
manuscript is to study how the ECS and related prop- 
erties such as the roughening temperature behave with 
increasing number of states. Like in two dimensions the 
model undergoes a phase transition at a temperature we 
will refer to as the melting temperature T m . Below this 
temperature the model has q different phases, each of 
which is dominated by one of the q possible spin values, 
whereas for temperatures above T m there is only a sin- 
gle phase in which on average all of the spin values are 
present in equal amounts. It is known that for q > 3 
the melting transition is of first order |l3|| , in contrast 
to the two-state (Ising) case where the transition is con- 
tinuous. It is therefore not clear a priori whether there 
is a roughening transition when q > 3, nor in how far 
the ECS should approach a sphere when T approaches 
Tr. Furthermore, one should expect the surface tension 
to approach a non-zero limit as the melting temperature 
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is approached from below, whereas for continuous tran- 
sitions it vanishes in this limit. 

The manuscript is organized as follows. First, in sec- 
tion |n| we outline the numerical procedure to efficiently 
determine the ECS and show the resulting shapes for all 
numbers of states for which accurate numerical estimates 
of the melting temperature have been reported: q = 2 
(Ising), 3,4 and 6. In section III we present some the- 
oretical considerations on the behavior of the model for 



large values of q. Section IV covers the measurements on 
the location of the roughening transition. In section |V] 
we present the orientation-dependent surface tension as 
extracted from the data on the ECS. We conclude with 
a discussion of the results and an outlook to further re- 
search. 



II. OBTAINING THE ECS 

Equilibrium shapes can be determined numerically for 
the g-states Potts model by forcing it into coexistence. 
This may be achieved by enforcing a constraint on the 
spin densities. The richest behavior is observed in case 
the maximum of q—1 constraints are enforced, but mostly 
the resulting configurations are still determined by the 
orientation dependence of the surface tension, and can 
thus be obtained indirectly from numerical simulations 
with a single constraint enforced. 

The constraint enforced in our simulations is the con- 
servation of the density of state 1, while the densities of 
all other states are allowed to fluctuate freely. To im- 
plement this, we combined a density-conserving cluster 
algorithm, described in a recent article |Q, with Glauber 
dynamics Q, constrained to never flip spins into or out 
of state 1. We chose the ratio of cluster updates and 
Glauber updates per site to be unity (i.e., a comparable 
amount of computational effort was spent to each) . The 
run-times of the simulations varied between 400 000 and 
800 000 Monte Carlo updates per site and we typically al- 
lowed the system to thermalize for 60 000 to 100 000 time 
steps. The lateral system size was 50 sites in all lattice 
directions, with periodic boundary conditions, and the 
fraction of conserved spins was typically 0.25. In order to 
directly obtain the ECS we repeatedly took "snapshots" 
of the evolving system, each time centering the cluster 
around the origin. Taking the autocorrelation time into 
account, we typically obtained 2000 independent equilib- 
rium snapshots from one simulation. At the end we cal- 
culated the time-averaged density profile of the conserved 
component. Our results for the ECS were obtained from 
the 50% iso-density surface in this time-averaged density 
profile. 

In figure ^, we show the equilibrium crystal shapes for 
several temperatures, for the case of q = 2 (Ising), 3, 
4, and 6 states, respectively. These shapes are obtained 
from the 50% iso-density surface in the density profile, 
after full symmetrization, i.e. after we averaged over all 



48 possible invariant mirror images of the cube. 

The simulation temperatures are chosen as fractions 
of the estimates for the melting temperatures reported 
in Refs. |T^-|T7|; these values are given in table ||. 



v« 4» 

v» v* v* 



3 




FIG. 1. Symmetrized Equilibrium Cluster Shapes. From 
top-left to bottom-right, the boxes show the shapes for the 2, 
3, 4 and 6-state Potts model, successively. Within each box, 
again from top- left to bottom-right, the temperature increases 
from 0.25 T m to 0.80 T m , in increments of 0.05 T m . 



III. BEHAVIOR FOR LARGE q 

For large values of q the thermal evolution of the ECS 
can be understood from the following arguments. First 
of all notice that at low temperatures the q phases of 
the q-state Potts model without conservation laws will 
be very close to the zero temperature phases in which all 
spins on the lattice have equal value. In d dimensions the 
free energy per site / of these states therefore satisfies 



! -dJ. 



(2) 
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The high temperature phase is completely dominated by 
the entropy resulting from the q different occupation op- 
tions for each site, with a resulting entropy per site of 
ks lng and free energy 

/ (Thigh) » -k B Tlnq. (3) 
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The melting temperature is obtained to a good approxi- 
mation by equating these two expressions, with the result 



2 lnq 



(4) 



for large q, with (3 — 1/fcgT. Note that for d = 2 this is 
in agreement with the exact result Jl^] (i m J — ln(l+y^). 
For large q one sees that J 3> ksT for temperatures be- 
low T m . Therefore local excitations, typically consisting 
of single overturned spins, will be extremely rare in the 
low temperature phases, as they require an energy far 
exceeding the gain in free energy due to the entropy in- 
crease. Similarly the tendency to form clusters of equal 
spin in the high temperature phase is very weak, as the 
resulting gain of energy by far does not compensate the 
resulting loss of entropic free energy. As a result the 
melting transition for high q will become very sharp. 

Now we can also estimate the orientation dependence 
of the surface tension. In good approximation, an inter- 
face of 001-orientation (or symmetrically equivalent) is 
just a flat interface between two pure states, with a free 
energy of J (one broken bond) per unit area. Similarly, 
by counting broken bonds, one finds that an interface of 
fcZl-orientation (here k and I are not necessarily integers) 
has an energy per unit area of (1 + k + O/Vl + fc 2 + I 2 . 
For non-zero temperature there is an additional contribu- 
tion from the entropy of the steps required to form such 
an interface. For not too large values of k and I this may 
be obtained by multiplying the step density needed to 
create a surface of orientation kll with the meandering 
entropy of such a step. The resulting expression is. 

+ I) J - k B T\n [(fc + l) k+l /(k k l 1 )] 
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This result is independent of the value of q, especially it 
also holds for the Ising model. Excitations creating spins 
not belonging to the two phases that coexist at the inter- 
face are so rare that they may be neglected completely. 
For k/l or l/k <C exp(— (3 J) the entropy becomes domi- 
nated by thermal step fluctuations, but this concerns an 
extremely small range of orientations only (However, it 
does set the distance between facet edges). Blote and 
Hilhorst jl8| give expressions for the free energy from 
which the low temperature ECS may be obtained for all 
orientations. 

The unimportance of all but the two coexisting phases 
implies that for large q the 001-facets do not roughen, 
cts Tjyi is well below the roughening temperature of these 
facets in the Ising model (see Section IV). The equilib- 
rium shape then remains nearly cubic for all tempera- 
tures up to T m . 



IV. LOCATION OF THE ROUGHENING 
TRANSITION 

Roughening transitions of crystal surfaces are char- 
acterized by the disappearance of macroscopically flat 



regions, or facets, in the ECS. Facets existing at low 
temperatures often disappear when the temperature in- 
creases, at a characteristic temperature known as the 
roughening temperature T R of the specific facet orien- 
tation. At this temperature the inverse radius of cur- 
vature (R c )~ l of the crystal surface at the center of 
the facet jumps from zero (when a facet is present) to 
a non-zero value. Using renormalization-group calcula- 
tions, Jayaprakash et al. Jl^,|l2| showed that the size of 
this jump satisfies the universal relationship 
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where zq is the distance from the tangent plane at the 
facet to the center of mass of the crystal, and 70 is the 
surface tension of the 001-facet at T — T R . The latter 
is not known exactly for any of the q-state Potts models 
in three dimensions. An approximation consisting of the 
ground state value plus the first correction term in a low 
temperature expansion is 



70 = (3 J - 2e 
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It will turn out that this approximation suffices for our 
goal: even for the smallest g-values the contribution of 
the correction term at T — T R is already smaller than 
the estimated systematic error in the curvature measure- 
ments. 

If the center of mass of the ECS is placed at the origin, 
the centers of the facets that are present below the rough- 
ening temperature are located on the principal axes. To 
estimate the location of the roughening transition, we 
first measure the curvature at the six points where the 
ECS (obtained as in the previous section) cuts a principal 
axis, in the two principal directions tangential to the sur- 
face. The curvature was obtained by fitting a quadratic 
function y = yo — ij/ix 2 ; the fitted result for yi/yo is 
the normalized inverse radius of curvature. Since we use 
a non-zero fitting range, the slope of the curve cannot 
become infinite and consequently we do not observe a 
jump in curvature; in fact we find quite smooth depen- 
dence of curvature on temperature, as a result of finite 
size, round-off and fluctuation effects. However, we made 
sure that our procedure to estimate Tr is not very sen- 
sitive to this. Notably our estimate of the roughening 
temperature for the Ising model is close to previous esti- 
mates: Adler reported T R /T, m = 0.55 ±0.02 [§0|, Mon et 
al. reported T R /T m = 0.54 ± 0.02 and Holzer and 
Wortis reported T R /T m = 0.545 ± 0.004 §§. 

The resulting measurements of the normalized inverse 
curvature are plotted in figure |. The error bars indi- 
cate the statistical error; additional simulations for some 
points indicate that the systematic errors arising from 
the effects of thermalization, cluster shape deformations, 
long lived thermal excitations, etc., are larger. In the 
same plot, the curves described by Eq. @ are also plot- 
ted. The intersection points of these two curves are esti- 
mates of the location of the roughening transition. The 
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resulting values for the roughening temperature can be 
found in table |; the error indications are our estimates 
of the statistical and systematic errors combined. 

In the previous section we noted that for large q there 
is no roughening transition. Obviously q — 6 is not large 
enough to observe this, but the increase of Tn/T m is very 
clear. It would be interesting to know the largest g-value 
for which a stable roughening transition does exist. The 
values of (3Jr quoted in table [| combined with the esti- 
mate of T m of Eq. (||) suggest a value between 20 and 
40. With the present method this would make simula- 
tions very slow, and it would also require an accurate 
measurement of the melting temperatures, since no liter- 
ature values exist. 

Roughening temperatures, g-state Potts-model 
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3 |ie] 

4 iyi 

6 



Tr/T„. 



0.84(2) 
0.86(2) 
0.89(2) 
0.93(2) 



0.443309(2) 
0.5505(1) 
0.631(2) 
0.751(2) 



0.53(3) 
0.64(3) 
0.71(3) 
0.81(3) 



0.8813 
1.0051 
1.0986 
1.2382 



TABLE I. Inverse transition temperature, ratio of the 
roughening and transition temperatures and inverse roughen- 
ing temperature of the q-states Potts model. The references 
in the first column point to the sources we used for the data 
on the transition temperatures. For comparison we added the 
critical interaction parameter j3J m = log(l + y 7 ?) for the two 
dimensional Potts model in the last column. 




0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 
T/T m 

FIG. 2. Normalized inverse curvature for the Potts model 
with q — 2 (pluses), 3 (crosses), 4 (stars) and 6 states 
(squares). The intersection points of the lines, described by 
Eq. (^|), and the data points are estimates for the roughening 
temperature Tr. 



V. ORIENTATION DEPENDENCE OF THE 
SURFACE TENSION 



In general, the surface tension 7 can be written as 7 = 
70 7 (n), where 70 is the surface tension of an interface 
oriented in the 001, or symmetrically equivalent direction 
]l2| . In experiments with Pb crystals in equilibrium with 
their vapor, Heyraud and Metois used the inverse Wulff 
construction (see Rcf. ) to obtain the angular part of 
the surface tension 7(77) as a function of orientation n and 
the same method has been employed by Surnev et al. [|l0| . 
Here we use it to obtain the orientation dependence of 
the surface tension for each ECS as obtained in section 

i 

The surface tension 7(77) is proportional to the dis- 
tance from the center of the cluster to the tangent plane 
perpendicular to n, touching the iso-surface in the point 
X. If the ECS is scaled such that the distance to the sur- 
face from the center of the shape along the lattice axes 
is unity, then 

7(77,) = max 7i ■ X, (8) 

where X is an element of the scaled iso-surface /. Thus, 
for every direction 77., we have to find the point X in the 
iso-surface, for which the number 77 • X is maximal. This 
procedure works in both two and three dimensions. 

Using ECS's obtained from fully symmetrized density 
profiles, we measured 7(77.) along the arc in the llfc-zone, 
which connects the 001, 111 and 110 directions. Figure || 
shows the results. 
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that aspheric deviations become larger with increasing q. 
As we saw in Section III the ECS becomes nearly cubic 
for large q. Obviously q — 6 still should be considered a 
small q- value in this context. 

If the ECS is not faceted, the angular part of the sur- 
face tension will approach the point 8 = with zero 
slope. If however the ECS has 001 facets, the approach 
to the point 8 = will occur with a non-zero slope, 
resulting in a cusp in the slope along the arc through 
this point [Q . Looking for the temperature where these 
cusps first appear is an alternative way to measure the 
roughening temperature. We found this to be less ac- 
curate than the procedure described in section IV. The 
results were however consistent. 



VI. DISCUSSION AND FUTURE RESEARCH 

We studied the Potts model with q = 2, 3, 4, and 6 
states, forced into coexistence by fixing the density of 
one state. The resulting ECS was studied as a function 
of temperature and number of states. 

We found that the roughening transition, which is well 
known for the Ising model (equivalent to the two-states 
Potts model), persists for higher numbers of states, at 
least up to six states. The temperature Tr at which this 
roughening transition occurs, measured as a fraction of 
the melting temperature, tends to increase with increas- 
ing number of states. 

In the future, we want to study the more general and 
richer behavior of the ECS in case more than one quantity 
is conserved. For instance, in the three-state Potts model 
close to its melting point, with the constraint p\ 3> P2 > 
ps, the ECS resembles the shape of two soap bubbles 
with a common interface; this changes under variations 
of temperature and the ratio pij p^. 

We also want to look at equilibrium shapes in con- 
strained geometries, like fluids between parallel plates, 
or systems with grain boundaries. 

Finally we are investigating the behavior for larger q- 
values with the aid of different Monte Carlo techniques. 



FIG. 3. Angular part of the surface tension 7(n) as a func- 
tion of orientation h, measured along the arc that connects 
the 001 (at 9 = 0), 111 (at 9 = 0.955) and 110 directions 
(at O = 7r/2, where 9 is the azimuthal angle). From top to 
bottom the plots show the 2, 3, 4 and 6 states results. In each 
plot, from the uppermost curve downwards, the temperature 
ranges from 0.25 T m to 0.80 T m , in increments of 0.05 T m . 

These figures show that the angular dependence of 
the surface tension becomes nearly constant at increas- 
ing temperature, for all numbers of states q considered 
here. Thus, for q < 6 the ECS approaches a sphere. For 
the Ising model (q = 2) this was to be expected, since 
this model undergoes a continuous phase transition at 
the melting point. On the other hand, it is also clear 
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